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Abstract 

We present a discrete stochastic model which represents many of 
the salient features of the biological process of wound healing. The 
model describes fronts of cells invading a wound. We have numerical 
results in one and two dimensions. In one dimension we can give 
analytic results for the front speed as a power series expansion in a 
parameter, p, that gives the relative size of proliferation and diffu- 
sion processes for the invading cells. In two dimensions the model 
becomes the Eden model for p ft; 1. In both one and two dimensions 
for small p, front propagation for this model should approach that 
of the Fisher-Kolmogorov equation. However, as in other cases, this 
discrete model approaches Fisher-Kolmogorov behavior slowly. 

1 Introduction 

The biology of wound healing is fairly well understood ^01 • A simplified 
version of the process may be given as follows: a layer of undamaged cells 
is usually quiescent, so that the birth rate of cells matches the death rate, 
and both are quite small. When a wound is suffered, there is a rapid signal 
the wakes the cells up - perhaps a pulse of ATP or a calcium wave. Cells 
at the edge of the wound become more mobile, and also enhance their 
proliferation rate. (Otherwise the healed layer would not have the right 
density.) A typical experiment to study this process consists in plating 
suitable cells (e.g. epithelial cells) on a substrate so that they form a 
confluent monolayer. Then a scratch is made in the layer, and the process 
of filling in the scratch is studied. For example, the speed of advance of 
the invading cells, v, is easily measured. 
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There have been many modeling studies of wound heahng [T^ini lSllTl ] . 

In many cases is an exception) the process is studied using some vari- 
ant of the Fisher-Kolmogorov (FK) equation 0] [7|. This is an obvious 
model to use. It builds in diffusion with diffusion constant D and prolifer- 
ation with growth rate k (related to inverse doubling time). It also shuts 
off growth for the confluent layer at density Co- 
de/ dt = DV^c + kc{l - c/co) (1) 

The justification for using a continuum equation for a cellular process relies 
on the common experience that coarse-graining is reasonable for dynamic 
processes involving a large number of agents. In this particular case, we 
expect that the FK equation should be useful if the characteristic length 
of the pattern predicted by Eq. 0J is much larger that the size of a cell. 

However, it is well known [210111 that coarse-graining the FK equation 
has many pitfalls even in this limit, and that the transition to the contin- 
uum limit is often very slow. This motivates the present investigation: we 
present a discrete stochastic model for wound healing, and study it in var- 
ious limiting regimes. It is quite similar to a model previously introduced 
and studied for flame- front propagation ^ISl. Thus, our results and meth- 
ods should be of interest beyond the explicit biological context. We will 
give new numerical and analytical results, and show how, in one and two 
dimensions, our model aproaches the FK limit. We will show that in the 
biologically relevant regime there are corrections to FK due to discreteness. 

2 Formulation of the model and known prop- 
erties 

Consider a set of sites that form a linear or square lattice, corresponding 
to one or two dimensional 'tissues'. We allow each site to be occupied by 
zero or one cells. Our initial configuration is an occupied half space: if 
i labels the x coordinates of the sites, then we have all sites with i < 
occupied. The dynamical rules are as follows: we choose a parameter p 
which specifies the proliferation rate of the cells. Then at any time step 
we choose a cell at random, and an adjacent site at random as a target for 
diffusion or proliferation. E.g., in Id if we choose a cell at site i, we also 
pick site i -|- 1 or i — 1 as a target . If the target site is empty, with probability 
p we put a new cell at the target, and with probability q = 1 — p we move 
the chosen cell to the target. If the target site is filled, we do nothing. 
These are examples of the elementary processes allowed: 

• (...1111000...) (...1111100...); probability p 
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• (...1111000...) ^ (...1110100..); probability g 

As time advances cells appear for i > 0. These form a front or chemical 
wave. We will examine the speed, v(j)), and front width, w{p) for the 
invading cells. Precise definitions for these quantities will be given below. 

An essentially identical model was devised by Kerstein [S] to describe 
flame-front propagation. He studied it numerically in Id, and Bramson et 
al. PP found some analytic results, also in Id. In their formulation there is 
a parameter 7 which may be identified as (1 — p)/p in our notation. Also, 
in their model the time unit is different from ours by a factor 1 + 7. If ^^(7) 
denotes the front speed in the Kerstein model, we have: 

v{p)=V{^)/{l + j). (2) 

In ^ there are two exact results. In our notation these are: 

• v{p) 1/2 + 0{q^) asp^l, 

• v{p) — > y/2p as p ^ 0. 

The first of these two is obvious. In the limit p — I there is no diffusion, 
only proliferation, and the half space advances with no vacancies. The only 
process allowed is to choose the leading cell at site i and proliferate at site 
i + 1. Since half the moves are wasted by choosing as a target the filled 
site at i — 1 the front speed is 1/2. The lack of a term linear in q will be 
derived below. 

The second result may be understood by comparison with Eq. Q . Con- 
sider the coarse-grained limit of our discrete model using the lattice con- 
stant as the unit of space, and a computer time step as the time unit. It 
is elementary to see that the diffusion coefficient, D, is 1/2. Now consider 
a collection of cells distant from one another with concentration c. In unit 
time the number will increase to {l+p)c. By integrating Eq. over space 
in the low density limit, we see that we must identify k — p. The front 
velocity given by Eq. ^ is well known for bounded initial conditions |11| : 

V = 2VDk = V2p- (3) 

Kerstein 5 verified both limits numerically. We will extend these one 
dimensional results below both numerically and analytically, and also in- 
vestigate the two-dimensional case. 

We note for future reference that the solutions of Eq. ^ generate an 
interface with an intrinsic width (see Figure 1^) given by: 

w = y^Djk cx 1/Vp. (4) 

Previous authors have not discussed the front width, but, as we will see, it 
is relevant to a biological interpretation of the results. 
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Figure 1: Sketch of traveling wave solution to the FK equation. The front 
position and width can be defined as shown. A is the negative derivative 
of c, see text. 



3 Numerical results 
3.1 Defining the front 

The solution to Eq. is a traveling front of the general form shown in 
the sketch in Figure ^ Our data for the discrete model is the form of 
occupancies of sites as a function of time. We present here a useful way to 
analyze such data that allows easy comparison to continuum theories. 

We start by defining the occupancy of a given column of our numerical 
data, P{i). In Id this is simply 1 or 0, depending on whether site i is 
occupied. In 2d it is the average occupancy of column i, that is, the number 
of occupied sites with first coordinate i divided by the width of the system 
(the total number of such sites) which we denote by L. We also define the 
negative of the discrete derivative of P; it is localized near the interface: 



Note that at long times we certainly have P{0) = 1, and for large enough 
i,P(i) = 0. Thus: 



A{i) ^ P{i) - P{i + 1). 



(5) 



OC 




(6) 
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That is, we can use A as a weight function to define averages. We put: 



2 = 1 = 1 

= f2^'Ai^)^Y^i2^-l)Pi^),■■■ (7) 



Here, Up is the number of particles for z > 1 in Id, or that number divided 
by L in 2d. That is, we get the position of the front by the total mass of 
created particles. 

The front speed is defined as 



V — lim (i) /t. 

t — >oo 



The front width is given by 



Other moments of the distribution can be defined similarly. 



3.2 One dimension 

The results of our simulations arc shown in Figures l(2Jl and (PJ . The front 
speeds were found by fitting (i) to vt. It is remarkable that the front speed 
is quite well defined even for very small p. Of course, for p = Q the speed 
is not defined at all. 

The convergence to the continuum predictions is quite evident in the 
figures. Note that the prediction for v{p) does not contain an adjustable 
constant, so the agreement is quite remarkable. However, following the 
work of [Sl|2j we would expect the corrections to the continuum prediction, 
Eq. ©, to follow: 

v{p)^^^-A/\u\p), (8) 

where A is a factor of order unity. In fact, this expression does not fit our 
results. Rather, the correction to the continuum formula is more like a 
power law with a power near 2/3. 



3.3 Two dimensions 

Our numerical results for v{p) in 2d are given in Figure^ Note that as 
p ^ 1, v{p) > 0.5. Analysis of the processes the contribute to front motion 
in 2d is more complex than in Id where v{l) = 0.5 is an exact result. 
In particular the front will always be rough (see below) so that particles 
behind the leading particle will not be blocked from advancing. 
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Figure 2: The front speed, v{p) in one dimension. The numerical simula- 
tions are averaged over 50 realizations for large p and up to 1000 for the 
smallest p to give the errorbars. Upper line is the continuum approxima- 
tion, Eq. Q. 
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Figure 3: The front width, w{p) from the same simulations as Fig. Q. 
The width is arbitrary up to a numerical factor. Hence the continuum 
approximation (upper line), Eq. Q is multiplied by a fitting factor. 
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Figure 4: The front speed in two dimensions. The continuum approxima- 
tion is from Eq. 



For p <C 1 we expect that we should converge to the resuh of the FK 
equation, namely that v{p) oc ^yp. However, for 2d we have no convincing a 
priori estimate of the pref actor. Following the treatment above, we might 
proceed by noting that in 2d = 1/4 for the discrete model. We then 
have: 

v{p) « VP- (9) 

As we can see from the figure, this is a reasonable estimate for small p. 

In two dimensions the width of the interface is a more complicated 
object than in one dimension . The reason for this is that in the presence 
of fluctuations the front can do two different things: it can spread so that 
it has an intrinsic width (as in Id) by having a reduced density in the 
interface region, but also it can wander. Indeed, for p = 1 wandering is the 
only effect possible. (Recall that in Id w{p = 1) = 0.) In fact, in the large 
p regime this model is identical to the Eden model |3] where perimeter sites 
all grow with equal probability. 

The phenomenology of the Eden model is well understood ^I]. The 
wandering of the interface is time-dependent and obeys (in our system of 
units): 



w cx t^/^ t < L^/^ 



(10) 




Figure 5: Scaling of w'^{p — 1) with time. From bottom to top, L = 
64, 128, 256, 512. Also shown are lines giving the expected scaling for early 
times, w'^ oc t'^^'^, and late times, w'^ oc L. 




Figure 6: The saturated value of w'^ in 2d. Also shown is the continuum 
approximation, w'^ cx 1/p. 
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This is indeed the case here: see Figure|31 The agreement with the Eq. pU|) 
is reasonable. We have verified that the scaling behavior given in Eq. H10() 
persists down to p = 0.5. 

However, as p decreases the intrinsic width grows rapidly. As soon 
as the intrinsic width exceeds the saturated width from wandering (the 
second line of Eq. CO))) we will loose the power-law time dependence of w. 
In Figure we show the saturated width for a range of p. 

4 Results for ^ « 1 in one dimension 

For p « 1 the dominant process is proliferation. For p = 1 this gives rise to a 
simple configuration as we have mentioned above: all sites behind the front 
are occupied, and the front advances because the leading cell proliferates. 
For q <^ 1 there is a small probability q/2 oi creating a configuration with 
a 'hole'. Because the model is very simple we can use this observation to 
work out the power series expansion of v{q). 

4.1 Exact solution of model for states with one hole 

Suppose we consider only states with zero holes or one hole at any position. 
We expect these to be the dominant configurations small q. Define the 
states: 

• |0) = (...11111000...) 

• |1) = (...11101000...) 

• |2) = (...11011000...) 

• |3) = (...10111000...), etc. 

We allow transitions only between these states. The transitions and their 
associated probabilities Wij = W{\i) are: 

Woo^p/2 Woi^q/2 

M^io-(l+p)/2 M^i2 = l/2 
W^o^p VK„,„_i=g/2 iy„,„+i = l/2 (n>l) (11) 

Note that in many of these transitions the actual location of the rightmost 
1 changes. We always define states in a frame moving with the front. 
The equations for the probabilities are: 

oo 
n=l 
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Pl(W^10 + Wx2) = PaWoi + P2W21 
-Pn(W^„o + Wn.n-1 + W„,„+l) = Pn-lWn-l,n + Pn+lWn+Un [u > 1) 

(12) 

Using Eq. (|ll|l we have: 

^'"'^Pn = (l)Pn-l + (i]P.H (n>l) (13) 



For n > 1, we make the ansatz P„ = a" ^Pi, and inserting this in the 
last equation above, we find: 



3-q-^9-10q + q^ 1 4q 2q^ 
" = 2q = 3 + 27 + 243---- ^^'^ 

Next, we substitute P2 = aPi into the second Une of Eq. (|13|l and use 
Sn°=2^n — 1 — -Pq — -Pi; in the first Une of Eq. (fT^ . Solving these two 
equations we find: 



Pn = 



Thus 



Pi 

P2 = aPi = 



2(l-g)(3-(l + a)g) 
6 - (5 + 2a)q + ag2 

^ q q^ n^ m^ 

2 12 216 3888 ■" ^ ' 



q q^ I9q'^ 



3 54 972 

9 ^ 162 ^ 2916 



The velocity can be found from 



l_^P(x01) 

2 2' ^ ^ 
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Figure 7: Results for v{p) in one dimension for 0.5 < p < 1. The o's are 
the numerical results, the dotted line is the quadratic approximation, and 
the solid line the power series of Eq. (|26(l . 



where x is any string of O's and I's, and we omit the zeros to the right. 
In this case, P(xOl) = Pi, because |1) is the only one-hole state that ends 
with (01). Thus, we have 

1 g2 _^ IV 

V — . . . (17) 

2 6 108 1944 ^ ' 

We have precise numerical data for v{q) for small q; see Figure [7| We 
find that Eq. ((T7|l is correct only up to quadratic order, as we might expect 
in the one-hole approximation. For example, for g = 0.1, we find numeri- 
cally that V = 0.498292, while 1/2 - = 0.49833, so that the coefficient 
of q'^ should be negative, not positive. Note q'^/108 « 0.00001. 

4.2 Reduced distribution functions. 

In order to go beyond quadratic terms in q, we introduce reduced distri- 
bution functions. This method would, in principle, allow the power-series 
expansion to be carried to arbitrary order. 

A reduced distribution function is probability to have a given pattern 
near the front for any pattern to the left. For example, as in Eq. (|16|l : 

P(xOl) = prob(. . . raxOlOO . . .) 
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where the sites marked as x are any string. Likewise, we define P(xll), 
P(xOOl), P(xlOl), and so forth. Note that, for example: 

P(xOOl) + P(xl01) = P(xOl) 
P(xOll) + P(xlll) = P(xll). 

We can derive a hierarchy of equations based on events that change the 
last n sites. For n — 2, consider all events that change the probability that 
the last two sites are (11). We have: 

' P P\ nl nni / 1 1 P 



l + JlflxOll) - {l\p{ym) = a. (IS) 



The positive terms represent events that increase the population of states 
ending with (11), and the negative terms represent events that decrease that 
population. Next we write all states in terms of P(xOl), P(xOOl), P(xOll), 
i.e., states that have a leading zero on the left. For the other states, we 
use: 

P(xlOl) ^ P(xOl) - P(xOOl) (19) 
P(xlll) = P(xll)-P(x011) = 1-P(x01)-P(x011). (20) 

Then Eq. lfTH|l becomes: 



2 V2 



)p(x01) + (^i±^^P(x001) + (I)P(XOII) = 0. (21) 



Now, we expect that P(xOl) ©(g), P(xOll) = 0(q), and P(xOOl) = 
C(g^), because a diffusive move (weight g/2) is required to produce each 
empty site starting from state |0). To order q we find from Eq. (|21|l 

P(x01) = | + O(a (22) 

which agrees with the leading behavior found in the one-hole approxima- 
tion. This implies that the velocity is given by: 

Note that there is no linear term, as mentioned above. 
For the next-order behavior, we use: 

P(xOll) - 1+0(9') (24) 
P(xOOl) = y+0(<z') (25) 



12 



where the leading behavior in Eq. H24|) is from the one-hole approximation. 
The second line follows from a simple argument: the leading behavior of 
P(xOOl) is determined by P(...1001), and its leading behavior is deter- 
mined by the equation: 

P(...1001)(3/2) = P(...101)(q/2) + P(...10001)(p/2 + q/2) + ... 

Again ... represents a string of all I's to the left. The second and higher- 
order terms on the right-hand-side are of order q'^ , so to leading order we 
find P(...1001) = (g/3)P(...101) = gV9, which proves Eq. 
Using these results, Eq. H21|l implies 

^^(x01) = | + ^+O(a 

which yields 

2 6 21 ' 
For the case q = 0.1, these three terms give v = 0.498296, in close agree- 
ment with the numerical simulations, which give 0.498292. 

We have carried this procedure to the next order, n = 3, by straightfor- 
ward extensions of what we have given above. The result for the velocity 
is: 

+ (26) 
2 6 27 1215 ^ > ^ > 

For q = 0.1, the predicted velocity is now 0.4982923, in complete agree- 
ment with the numerical result 0.498292. Even a,t q = 0.5, the prediction 
of Eq. (PI), w = 0.4511831, is within 0.2% of the measured value, 0.45014. 



5 Application to biology 

Our emphasis in this paper has been an analysis of the model introduced 
in the introduction. It is interesting, nevertheless, to make some comments 
on the relationship of this model to real biological systems. Needless to say, 
our view of wound repair is very much oversimplified. In a real tissue there 
are various types of cells such as stem cells which have different behavior 
with respect to proliferation than others. Further, the proliferation cycle is 
complex, and involves time delays that we have not considered except in a 
rough way. Also, the initiation of wound healing is probably mediated by 
chemical signals rather than cell proximity as we have assumed |14| . 

However, if we are interested in macroscopic features such as the ve- 
locity and shape of the moving front, we are entitled to hope that many 
of these details will be unimportant. We can then ask how to translate 
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the parameters of our model to a real system. We will take as an example 
the experiment of Sheardown and Cheng |12| on the wounding of rabbit 
corneas. 

In the emphasis was on modeling with the FK equation. To this 
end the authors measured D in Eq. ^ by looking at the initial stage of 
invasion of cultured cells, and found D — 1.61 • 10~^mm^/s. The parameter 
k in Eq. is related to the mitotic rate of cells which was measured by 
labehng with a dye: = 4.3 days. Using these parameters the authors 
found reasonable agreement for the velocity of the front. Further, for these 
parameters, the shape of the front is quite 'fuzzy', that is, the width, w, 
is many cells across so that the wound fills in gradually, as observed. We 
should note that this is in sharp contrast to other observations |13[I5] where 
the advancing front is quite sharp. We will return to this point below. 

We now attempt to translate these observations into the parameters of 
our discrete model. We need to define units of length and time. For length 
it is natural to take a typical cell size, d = 10/i as the lattice unit. It is 
clear that we can define the hopping time, Thop by D = (P/2Thop- This 
turns out to be about 108 seconds for the rabbit cornea. However, there 
is another characteristic time, the cell cycling time, Tcyc = l/k. This is 
3.7 • 10^ seconds for the same experiment. In our model, in TV computer 
cycles there are Np cell cycles and Nq hops. Thus our time unit should be 

To determine the biological p note that p/q = Tcyc/Thop- For the rabbit 
experiment we get p — i ■ lO^**. This is the regime of very diffuse, fuzzy 
interfaces, as observed. In this regime FK modeling should be reasonable, 
though, as Figure El shows, there are still differences between FK and the 
discrete model in this regime. 

If we apply the same set of considerations to the systems studied in 
| 13l H] we find a contradiction. The width of the interface should be quite 
substantial for the small p's relevant to biological experiments, cf. Figure|31 
In fact, FK modeling shows the same thing. However, direct observation 
in these cases shows that the front is quite sharp. 

A possible solution to this quandary is given by jJJ where it is pointed 
out that cell-cell adhesion can have an effect on wound healing in some sys- 
tems. In fact, for the cells that they study they can regulate the adhesion, 
and hence the front width, by controlling the supply of Ca++ ion in the 
solution bathing the cells. A detailed discrete model in their paper shows 
these effects too. This is an interesting avenue for future work. 
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6 Summary and conclusions 



In this paper we have extended the work of on a discrete model. We 
have shown that the model can be interpreted as a representation of the 
important biological process of wound healing. We have given numerical 
results in one and two dimensions, and a power-series expansion of the ve- 
locity around q — in one dimension. We have shown that the biologically 
interesting regime is that of p ^ 1. 

There are a number of further extensions of this work that could be 
pursued. Our method of reduced distribution functions should be appli- 
cable to models with more complex rules as long as a sensible expansion 
parameter, analogous to q, is present. We do not understand why the con- 
vergence to the FK limit is different in our case than in the generic cases 
discussed in jUj. A more extensive numerical study may be called for. 

We think that the most interesting extension of the model would be to 
include cell-cell adhesion, in the spirit of ^1^]. This work is in progress. 
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